##### name of package:wkptest######

#### purpose: the Painleve test of polynomial PDE(s)####

##### date:2003-09-03


#####author: Gui-qiong Xu and Zhi-bin Li

wkptest:=module()
 export pltest;
 local  normeqs,diffrule,pbranchk,algbranchk,finda,kesi,
 extract_functions,initorder,seperatenl,leadorder1,
 leadorder2,leadorder3, leadanalysis, trunform,findres,
 simpltm0,createqs,checkcond,reducecond;
 option package;

 pltest:=proc(eqlist::list,funcs_list::list)
 global `wkptest/xi`,`wkptest/eta`,`wkptest/eq`,`wkptest/m`,
 `wkptest/lf3`,`wkptest/lx1`,`wkptest/lx2`,`wkptest/lx3`,
 `wkptest/lx4`,`wkptest/ltt`,`wkptest/ly`,`wkptest/lty`,
 `wkptest/eqnum`,`wkptest/settime`,`wkptest/nam`,`wkptest/paranum`,
 `wkptest/argf`,`wkptest/func`,`wkptest/passflag`;
 local cpu_time,i,ll1;
 `wkptest/settime`:=time();
 `wkptest/paranum`:=nargs;
 if `wkptest/paranum`=2 then
    `wkptest/argf`:=args[2];
 fi;
 `wkptest/m`:='`wkptest/m`';
 `wkptest/passflag`:=1;
 if nargs<1 then
   ERROR(` Too few arguments ! `);
 elif nargs>2 then
   ERROR(` Too many arguments ! `);
 elif nops(eqlist)<1 then
   ERROR(` You must specify a list of equations ! `);
 fi;
 for i from 1 to nops(eqlist) do
   if evalb(not(type(op(i,eqlist),equation))) then
     ERROR(` You should input equation forms!`);
   fi;
 od;
 `wkptest/eq`:=eqlist;
 `wkptest/eqnum`:=nops(`wkptest/eq`);
 printf(`The input equations are:`);
 for i from 1 to `wkptest/eqnum` do
   print(op(i,`wkptest/eq`));
 od;
 normeqs();
 leadanalysis();
 if `wkptest/passflag`=0 then
    ##2003-09-03
    ##printf(`The equation fails the Painleve test in leading analysis\n`);
    printf(`The leading orders and coefficients cannot be determined\n`);
    cpu_time:=time()-`wkptest/settime`;
    printf(`The running time is %a\n`,cpu_time);
    return;
 else
    trunform();
    findres();
    if `wkptest/passflag`=0 then
      ##2003-09-03
      ##printf(`The equation fails the Painleve test!\n`);
      cpu_time:=time()-`wkptest/settime`;
      printf(`The running time is %a\n`,cpu_time);
      return;
    else
      printf(`Please wait for a moment!\n`);
      createqs();
      checkcond();
      ##if `wkptest/passflag`=0 then
        ##printf(`The equation fails the Painleve test in checking conditions\n`);
        cpu_time:=time()-`wkptest/settime`;
        printf(`The running time is %a\n`,cpu_time);
        return;
      ##else
      ## printf(`The equation(s) passes(pass) the Painleve test\n`);
      ## cpu_time:=time()-`wkptest/settime`;
      ## printf(`The running time is %a\n`,cpu_time);
     ##fi;
    fi;
 fi;
 end:

 extract_functions:=proc()
 global `wkptest/lv1`,`wkptest/func`;
 local  ii,tset1,tset2,lw,temp;
 lw:=1;
 tset1:=`wkptest/lv1`;
 tset2:={NULL};
 while evalb(lw>0) do
   lw:=0;
   for ii from 1 to nops(tset1) do
     if type(op(ii,tset1),`+`) then
         tset2:=tset2 union {op(op(ii,tset1))};
         lw:=1;
     elif type(op(ii,tset1),`*`) then
         tset2:=tset2 union {op(op(ii,tset1))};
         lw:=1;
     elif type(op(ii,tset1),`^`) then
         tset2:=tset2 union {op(op(ii,tset1))};
         lw:=1;
     elif has(op(ii,tset1),diff) then
         tset2:=tset2 union {op(op(ii,tset1))};
         lw:=1;
     else
         tset2:=tset2 union {op(ii,tset1)};
     fi;
   od;
   tset1:=remove(type,tset2,constant);
   tset2:=remove(type,tset1,name);
   tset1:=tset2;
   tset2:={NULL};
 od;
 `wkptest/func`:=tset1;
 end:

 kesi:=proc(ll::list)
   global `wkptest/eta`;
   local ll1,ll2,j,i,lk1,ll3,k;
   k:=1;
   ll1:=[k,-k*l,-k*n,-k*h];
   ll3:=[];
   for j from 1 to nops(ll)-1 do
     ll3:=[op(ll3),op(j,ll1)];
   od;
   ll3:=[op(ll3),-k*c];
   ll2:=zip((x,y)->x*y,ll3,ll);
   `wkptest/eta`:=add(j,j=ll2);
 end:

 diffrule:=proc(indep::list)
 local ll0,ll1,ll2,ll3,ll4,ll5,ii;
 ll0:=indep;
 ll1:=op(1,indep);
 ll2:=remove(has,indep,ll1);
 ll3:=[diff(phi(op(ll0)),ll1)=1];
 for ii from 1 to nops(ll2) do
   ll4:=op(ii,ll2);
   ll5:=diff(phi(op(ll0)),ll4)=-diff(psi(op(ll2)),ll4);
   ll3:=[op(ll3),ll5];
 od;
 return ll3;
 end:

 pbranchk:=proc(resl::list)
 local ii,ll0,ll1,ll2,counter;
 ll0:=resl;
 counter:=0;
 ll1:=1;
 for ii from 1 to nops(ll0) do
   ll2:=op(ii,ll0);
   if type(ll2,nonnegint) then
      counter:=counter+1;
   fi;
 od;
 if counter<>(nops(ll0)-1) then
   ll1:=0;
 fi;
 return ll1;
 end:

 algbranchk:=proc(resl::list)
 local ii,ll0,ll1;
 ll0:=resl;
 ll1:=0;
 for ii from 1 to nops(ll0) do
   if type(op(ii,ll0),integer) then
      next;
   else
      ll1:=1;
      break;
   fi;
 od;
 return ll1;
 end:

normeqs:=proc()
 global `wkptest/eq`,`wkptest/eqnum`,`wkptest/func`,
 `wkptest/var1`,`wkptest/eta`,`wkptest/lv1`,`wkptest/lv`,
 `wkptest/argf`,`wkptest/paranum`;
 local ii,jj,kk,tflag,ll0,ll1,ll2,ll3,ll4,tm1,tm2,ttm1;
 unassign('`wkptest/lv`','`wkptest/lv1`');
 ll0:=`wkptest/eqnum`;
 ll1:=`wkptest/eq`;
 ttm1:=array(1..ll0);
 `wkptest/lv`:=array(1..nops(ll1));
 if `wkptest/paranum`=2 then
   `wkptest/func`:=`wkptest/argf`;
   `wkptest/var1`:=[op(op(1,`wkptest/func`))];
 for ii from 1 to ll0 do
   tm1:=[];
   tm2:=[];
   ll2:=expand(lhs(op(ii,ll1))-rhs(op(ii,ll1)));
   ttm1[ii]:=ll2;
   for jj from 1 to nops(ll2) do
     tflag:=0;
     ll3:=op(jj,ll2);
     ll4:=indets(ll3);
     for kk from 1 to ll0 do
       if member(op(kk,`wkptest/func`),ll4) then
         tflag:=1;
         break;
       fi;
     od;
     if tflag=1 then
       tm1:=[op(tm1),ll3];
     else
       tm2:=[op(tm2),ll3];
     fi;
    od;
    `wkptest/lv`[ii]:=add(pp,pp=tm1);
  od:
 else
   for ii from 1 to ll0 do
   ll2:=expand(lhs(op(ii,ll1))-rhs(op(ii,ll1)));
   ttm1[ii]:=ll2;
   ll3:=remove(has,ll2,name);
   ll4:=remove(has,ll3,constant);
   `wkptest/lv`[ii]:=ll4;
   od;
   `wkptest/lv1`:=convert(`wkptest/lv`[1],set);
   for ii from 1 to nops(ll1) do
     `wkptest/lv1`:=`wkptest/lv1` union convert(`wkptest/lv`[ii],set);
   od;
   extract_functions();
   `wkptest/var1`:=[op(op(1,`wkptest/func`))];
 fi;
 `wkptest/eq`:=ttm1;
 kesi(`wkptest/var1`):
 printf(`\n`);
 end:

initorder:=proc()
 global `wkptest/eqm`,`wkptest/eta`,`wkptest/xi`,`wkptest/ly`,
 `wkptest/lv`,`wkptest/m`,`wkptest/func`,`wkptest/funname`,
 `wkptest/eqnum`;
 local equ1,equ2,equ3,equ,lf1,lf2,lf3,ee,mx,we,lk,we2,
       i,j,ij,jj,mm,ss,ll0,name1;
 unassign('eqm','mm','ss');
 ll0:=`wkptest/eqnum`;
 for i from 1 to ll0 do
   lk:=[];
   if nops({op(`wkptest/func`[1])})>1 then
     for j from 1 to nops(`wkptest/lv`[i]) do
        equ:=op(j,`wkptest/lv`[i]);
        equ1:=eval(subs({seq(`wkptest/func`[ij]=op(0,
        `wkptest/func`[ij])(`wkptest/eta`)
             ,ij=1..ll0)},equ));
        equ2:=simplify(subs(`wkptest/eta`=`wkptest/xi`,equ1));
        equ3:=convert(equ2,diff,`wkptest/xi`);
        lk:=[op(lk),equ3];
     od;
   else
     lk:=`wkptest/lv`[i];
   fi;
   `wkptest/eqm`[i]:=add(jj,jj=lk);
 od;
 `wkptest/funname`:=[];
 for i from 1 to ll0 do
   mm[i]:=10^i;
   name1:=op(0,`wkptest/func`[i]);
   `wkptest/funname`:=[op(`wkptest/funname`),name1];
   ss[i]:=name1(`wkptest/xi`)=tan(`wkptest/xi`)^mm[i];
 od;
 `wkptest/ly`:=[];
 `wkptest/m`:='`wkptest/m`';
 for i from 1 to ll0 do
   lf1:=[op(`wkptest/eqm`[i])];
   lf2:=[];
   for j from 1 to nops(lf1) do
     ee:=expand(simplify(subs({seq(ss[jj],jj=1..ll0)},op(j,lf1))));
     if(nops([ee])=1) then mx:=degree(ee,tan(`wkptest/xi`));
     else
         mx:=max(seq(degree(j,tan(`wkptest/xi`)),j=ee));
     fi;
     lf2:=[mx,op(lf2)];
   od;
   lf3:=[];
   for j from 1 to nops(lf2) do
     we:=0;
     we2:=op(j,lf2);
     for jj from ll0 to  1 by -1 do
        we:=we+iquo(we2,mm[jj])*`wkptest/m`[jj];
        we2:=we2-iquo(we2,mm[jj])*mm[jj];
     od;
     we:=we+we2;
     lf3:=[we,op(lf3)];
   od;
   `wkptest/ly`:=[op(`wkptest/ly`),lf3];
 od:
 end:

leadorder1:=proc()
 global `wkptest/ly`,`wkptest/m`,`wkptest/lty`;
 local i,ii,hh,m1_max,m2_max,m3_max,m_ineq,m_eq,max_m_value,ll1,ll2;
 m1_max:=max(op(op(1,`wkptest/ly`)));
 m2_max:={convert(m1_max,piecewise)};
 m3_max:={op(op(m2_max))};
 m_ineq:={NULL};
 for i from 1 to nops(m3_max) do
   if not (type(op(i,m3_max),`+`) or type(op(i,m3_max),`*`)) then
      m_ineq:=m_ineq union {op(i,m3_max)};
   fi;
 od;
 m_ineq:=remove(type,m_ineq,constant);
 m_ineq:=remove(type,m_ineq,name);
 m_eq:=map(convert,m_ineq, equality);
 `wkptest/lty`:={NULL};
 for i from 1 to nops(m_eq) do
    if type(lhs(op(i,m_eq)),constant) then
      hh:=rhs(op(i,m_eq))=lhs(op(i,m_eq));
      `wkptest/lty`:=`wkptest/lty` union {hh};
    else
      `wkptest/lty`:=`wkptest/lty` union {op(i,m_eq)};
    fi;
 od;
 ll1:=select(type,`wkptest/lty`,equation);
 ll2:={};
 for ii from 1 to nops(ll1) do
   ll2:={op(ll2),{op(ii,ll1)}};
 od;
 `wkptest/lty`:=ll2;
 finda():
 end:

 seperatenl:=proc(lz::list)
 global `wkptest/eqnum`,`wkptest/lp`,`wkptest/lq`,`wkptest/m`;
 local i,j,jj,ltt,sum1,maxlist,lp11,lq11,ll0;
 unassign('lp11','lq11');
 ll0:=`wkptest/eqnum`;
 for i from 1 to ll0 do
    lp11[i]:=[];
    lq11[i]:=[];
 od;
 for i from 1 to ll0 do
    ltt:=convert(op(i,lz),set);
    for j from 1 to nops(ltt) do
       sum1:=0;
       for jj from 1 to ll0 do
             sum1:=sum1+quo(op(j,ltt),`wkptest/m`[jj],`wkptest/m`[jj]);
       od;
       if sum1=1 then
          lp11[i]:=[op(lp11[i]),op(j,ltt)];
       else
          lq11[i]:=[op(lq11[i]),op(j,ltt)];
       fi;
    od;
 od;
 for i from 1 to ll0 do
   if lp11[i]<>[] then
      maxlist[1]:=max(op(lp11[i]));
      if hasfun(maxlist[1],max) then
         lp11[i]:=[op(maxlist[1])];
      else
         lp11[i]:=[maxlist[1]];
      fi;
   else
      next;
   fi;
 od;
 for i from 1 to ll0 do
   if lq11[i]<>[] then
      maxlist[2]:=max(op(lq11[i]));
      if hasfun(maxlist[2],max) then
         lq11[i]:=[op(maxlist[2])];
      else
         lq11[i]:=[maxlist[2]];
      fi;
   else
      next;
   fi;
 od;
 for i from 1 to ll0 do
    `wkptest/lp`[i]:=lp11[i];
    `wkptest/lq`[i]:=lq11[i];
 od;
 end:

 leadorder2:=proc()
 global `wkptest/lp`,`wkptest/lq`,`wkptest/m`,`wkptest/lty`;
 local lc,ld,i,j,jj,ii,max_rem1,max_rem2,belist,lists,ltt,lsm,lts,lty1,
       sd,sv,sv1,sv2,lth,ll0,lp11,lq11,xty;
 ll0:=`wkptest/eqnum`;
 unassign('lp11','lq11');
 for i from 1 to ll0 do
   lp11[i]:=`wkptest/lp`[i];
   lq11[i]:=`wkptest/lq`[i];
 od;
 for i from 1 to ll0 do
   belist[i]:={NULL};
   if lp11[i]=[] then
      for j from 1 to nops(lq11[i])-1 do
         for jj from 2 to nops(lq11[i]) do
           if j<>jj then
              belist[i]:=belist[i] union {op(j,lq11[i])=op(jj,lq11[i])};
           fi;
         od;
      od;
   elif lq11[i]=[] then
      for j from 1 to nops(lp11[i])-1 do
         for jj from 2 to nops(lp11[i]) do
           if j<>jj then
              belist[i]:=belist[i] union {op(j,lp11[i])=op(jj,lp11[i])};
           fi;
         od;
      od;
   else
      for j from 1 to nops(lp11[i]) do
         for jj from 1 to nops(lq11[i]) do
           belist[i]:=belist[i] union {op(j,lp11[i])=op(jj,lq11[i])};
         od;
      od;
   fi;
 od;
 for i from 1 to ll0 do
   if nops(lq11[i])>1 then
     for j from 1 to nops(lq11[i])-1 do
         for jj from 2 to nops(lq11[i]) do
           if j<>jj then
             belist[i]:=belist[i] union {op(j,lq11[i])=op(jj,lq11[i])};
           fi;
         od;
      od;
   fi;
 od;
 lists[2]:={NULL};
 for i from 1 to nops(belist[1]) do
   for j from 1 to nops(belist[2]) do
      lists[2]:=lists[2] union {{op(i,belist[1]),op(j,belist[2])}};
   od;
 od;
 for i from 3 to ll0 do
    lists[i]:={NULL};
    for j from 1 to nops(lists[i-1]) do
       for jj from 1 to nops(belist[i]) do
         lists[i]:=lists[i] union
         {{op(op(j,lists[i-1])),op(jj,belist[i])}};
       od;
    od;
 od;
 lists[1]:={};
 for i from 1 to nops(lists[ll0]) do
    if hastype(op(i,lists[ll0]),`=`) then
         lists[1]:=lists[1] union op(i,lists[ll0]);
    fi;
 od;
 xty:=[];
 lsm:={seq(`wkptest/m`[j],j=1..ll0)};
 for i from 1 to nops(lists[ll0]) do
    ltt:=solve(op(i,lists[ll0]),lsm);
    xty:=[op(xty),ltt];
 od;
 xty:=convert(convert(xty,set),list);
 ltt:=[];
 for i from 1 to nops(xty) do
   if nops(op(i,xty))=ll0 then
     ltt:=[op(ltt),op(i,xty)];
   fi;
 od;
 xty:=ltt;
 ltt:={NULL};
 for i from 1 to nops(xty) do
   if hastype(op(i,xty),name=0) then
      ltt:=ltt union {op(i,xty)};
   fi;
 od;
 xty:=convert(xty,set) minus ltt;
 xty:=convert(xty,list);
 ltt:=[];
 for i from 1 to nops(xty) do
  lts:=select(hastype,op(i,xty),{name=numeric,name=name});
  ltt:=[op(ltt),lts];
 od;
 xty:=ltt;
 lty1:=[];
 for i from 1 to nops(xty) do
  lts:=select(hastype,op(i,xty),name=numeric);
  lty1:=[op(lty1),lts];
 od;
 ltt:=[];
 for i from 1 to nops(xty) do
   if nops(op(i,xty))=ll0 then
     ltt:=[op(ltt),op(i,xty)];
   fi;
 od;
 xty:=ltt;
 ltt:=[];
 for i from 1 to nops(lty1) do
   if nops(op(i,lty1))=ll0 then
      ltt:=[op(ltt),op(i,lty1)];
   fi;
 od;
 lty1:=ltt;
 xty:=convert(xty,set);
 lty1:=convert(lty1,set);
 xty:=xty minus lty1;
 lth:={};
 for i from 1 to nops(xty) do
   ltt:=select(hastype,op(i,xty),name=name);
   if nops(ltt)>1 then
       ltt:={op(1,ltt)};
   fi;
   sv1:={};
   for j from 1 to ll0 do
       lp11[j]:=simplify(subs(op(i,xty),lp11[j]));
       lq11[j]:=simplify(subs(op(i,xty),lq11[i]));
       sv[j]:={solve(max(op(lp11[j]))=max(op(lq11[j])),{lhs(op(ltt))})};
       if evalb(sv[j]<>{NULL} and hastype(sv[j],
          {name=numeric,name<=numeric})) then
          sv1:=sv1 union sv[j];
       fi;
   od;
   sv2:={};
   for j from 1 to nops(sv1) do
      if rhs(op(op(j,sv1)))=0 then
          sv2:=sv2 union {op(j,sv1)};
      fi;
   od;
   sv1:=sv1 minus sv2;
   sv2:={};
   for j from 1 to nops(sv1) do
      if type(op(j,sv1),`=`) then
         lth:=lth union {{op(select(hastype,op(i,xty),
         name=numeric)), op(j,sv1)}};
     else
         for jj from 1 to rhs(op(op(j,sv1))) do
            lth:=lth union {{op(select(hastype,op(i,xty),
            name=numeric)), lhs(op(op(j,sv1)))=jj}} ;
         od;
     fi;
   od;
 od;
 xty:=lty1 union lth;
 lth:={};
 for i from 1 to nops(xty) do
    if hastype(op(i,xty),name=0) then
      lth:=lth union {op(i,xty)};
    fi;
 od;
 xty:=xty minus lth;
 lth:=[];
 for i from 1 to nops(xty) do
   if nops(op(i,xty))=ll0 then
      lth:=[op(lth),op(i,xty)];
   fi;
 od;
 xty:=lth;
 `wkptest/lty`:=xty;
 finda():
 end:

leadorder3:=proc()
 global  `wkptest/eqnum`, `wkptest/lp`, `wkptest/lq`, `wkptest/m`, `wkptest/lty`;
 local
 lc,ld,i,j,jj,ii,max_rem1,max_rem2,belist,lists,ltt,lsm,lts,lty1,
       sd,sv,sv1,sv2,lth,ij,ji,lvar,lvar1,lis1,lw,luy,luy1,ll0,xty,lp11,lq11;
 unassign('lp11','lq11','xty');
ll0:= `wkptest/eqnum`;
 for i from 1 to ll0 do
   lp11[i]:= `wkptest/lp`[i];
   lq11[i]:= `wkptest/lq`[i];
 od;
  for i from 1 to ll0 do
   belist[i]:={NULL};
   if lp11[i]=[] then
      for j from 1 to nops(lq11[i])-1 do
         for jj from 2 to nops(lq11[i]) do
           if j<>jj then
              belist[i]:=belist[i] union {op(j,lq11[i])=op(jj,lq11[i])};
           fi;
         od;
      od;
   elif lq11[i]=[] then
      for j from 1 to nops(lp11[i])-1 do
         for jj from 2 to nops(lp11[i]) do
           if j<>jj then
              belist[i]:=belist[i] union {op(j,lp11[i])=op(jj,lp11[i])};
           fi;
         od;
      od;
   else
      for j from 1 to nops(lp11[i]) do
         for jj from 1 to nops(lq11[i]) do
           belist[i]:=belist[i] union {op(j,lp11[i])=op(jj,lq11[i])};
         od;
      od;
   fi;
 od;
 for i from 1 to ll0 do
   if nops(lq11[i])>1 then
     for j from 1 to nops(lq11[i])-1 do
         for jj from 2 to nops(lq11[i]) do
           if j<>jj then
             belist[i]:=belist[i] union {op(j,lq11[i])=op(jj,lq11[i])};
           fi;
         od;
      od;
   fi;
 od;
 lists[2]:={NULL};
 for i from 1 to nops(belist[1]) do
   for j from 1 to nops(belist[2]) do
      lists[2]:=lists[2] union {{op(i,belist[1]),op(j,belist[2])}};
   od;
 od;
 for i from 3 to ll0 do
    lists[i]:={NULL};
    for j from 1 to nops(lists[i-1]) do
       for jj from 1 to nops(belist[i]) do
         lists[i]:=lists[i] union
         {{op(op(j,lists[i-1])),op(jj,belist[i])}};
       od;
    od;
 od;
 xty:=[];
 for i from 1 to 3 do
   for j from  1 to 3 do
      for ii from 1 to 3 do
        for jj from 1 to 3 do
           lvar:={ `wkptest/m`[1]=i, `wkptest/m`[2]=j, `wkptest/m`[3]=ii, `wkptest/m`[4]=jj};
           lvar1:=lvar;
           for ij from 1 to nops(lvar) do
              if evalb(op(lhs(op(ij,lvar)))>ll0) then
                  lvar1:=lvar1 minus {op(ij,lvar)};
              fi;
           od;
           lvar:=lvar1;
           lis1:=simplify(subs(lvar,lists[ll0]));
           for ij from 1 to nops(lis1) do
              lw:=0;
              luy:=op(ij,lis1);
              for ji from 1 to nops(luy) do
                 if not(evalb(op(ji,luy)))  then
                    lw:=1;
                 fi;
              od;
              if lw=0 then
                 xty:=[op(xty),lvar];
                 break;
              fi;
           od;
        od;
      od;
   od;
 od;
 lth:=convert(xty,set);
 xty:=lth;
 for i from 1 to nops(xty) do
   luy:=op(i,xty);
   luy1:=[];
   for j from 1 to nops(luy) do
     luy1:=[op(luy1),rhs(op(j,luy))];
   od;
   if numboccur(luy1,3)>=2 then
       lth:=lth minus {luy};
   fi;
 od;
 xty:=convert(lth,list);
 `wkptest/lty`:=xty;
 finda():
 end:

 finda:=proc()
 global `wkptest/eqnum`,`wkptest/m`,`wkptest/lty`,`wkptest/ly`,
 `wkptest/lord`,`wkptest/lpw`;
 local ii,jj,tm1,tm2,ll0,ll1,ll2,ll3,lr1,lr2,lr3,lr4;
 tm1:=[];
 tm2:=[];
 ll0:=`wkptest/eqnum`;
 ll1:=`wkptest/lty`;
 ll2:=`wkptest/ly`;
 for ii from 1 to nops(ll1) do
   ll3:=op(ii,ll1);
   lr1:=[];
   for jj from 1 to ll0 do
     lr2:=select(has,ll3,`wkptest/m`[jj]);
     lr3:=-rhs(op(lr2));
     lr4:=alpha[jj]=lr3;
     lr1:=[op(lr1),lr4];
   od;
   lr2:=[];
   for jj from 1 to ll0 do
     lr3:=op(jj,`wkptest/ly`);
     lr4:=-max(op(subs(ll3,lr3)));
     lr2:=[op(lr2),lr4];
   od;
   tm1:=[op(tm1),lr1];
   tm2:=[op(tm2),lr2];
 od;
 `wkptest/lord`:=tm1;
 `wkptest/lpw`:=tm2:
 printf(`\n`);
 end:

leadanalysis:=proc()
 global `wkptest/eqnum`,`wkptest/var1`,`wkptest/func`,`wkptest/lord`,
 `wkptest/lpw`,`wkptest/ly`,`wkptest/lty`,`wkptest/ltm0`,`wkptest/passflag`,`wkptest/lv`;
 local ii,jj,kk,rulef,eqcoefs,ll0,ll1,ll2,ll3,ll4,ll5,ll6,lr0,lr1,lr2,ls0,ls1,ls2,tm0,tm1,tm2,lpflag,lp1;
 ll0:=`wkptest/eqnum`;
 lr0:=`wkptest/var1`;
 ls0:=`wkptest/func`;
 initorder();
 if `wkptest/eqnum`=1 then
    leadorder1();
 elif `wkptest/eqnum`=2 then
    seperatenl(`wkptest/ly`);
    leadorder2();
 elif `wkptest/eqnum`>2 then
    seperatenl(`wkptest/ly`);
    leadorder3();
 fi;
 lr1:=`wkptest/lord`;
 lr2:=`wkptest/lpw`;
 tm1:=[];
 tm2:=[];
 tm0:=[];
 if nops(lr1)=0 then
   `wkptest/passflag`:=0;
   return;
 fi;
 for ii from 1 to nops(lr1) do
   ll1:=op(ii,lr1);
   #for the sake all orders should be negative
   lpflag:=1;
   for jj from 1 to ll0 do
     lp1:=rhs(op(jj,ll1));
     if not (type(lp1,negint)) then
        lpflag:=0;
        break;
     fi;
   od;
   if lpflag=0 then
      next;
   fi;
   rulef:=[];
   for jj from 1 to ll0 do
     ll2:=rhs(op(jj,ll1));
     ll3:=op(jj,ls0);
     ll4:=ll3=op(0,ll3)[0](op(lr0))*phi(op(lr0))^(ll2);
     rulef:=[op(rulef),ll4];
   od;
   eqcoefs:=[];
   for jj from 1 to ll0 do
     ll2:=expand(`wkptest/lv`[jj]);
     ll3:=op(jj,op(ii,lr2));
     ll4:=expand(eval(subs(rulef,ll2)*phi(op(lr0))^(-ll3+1)));
     ll5:=combine(ll4,power);
     ls1:=[];
     for kk from 1 to nops(ll5) do
       ls2:=op(kk,ll5);
       if member(phi(op(lr0)),{op(ls2)}) then
          ls1:=[op(ls1),ls2];
       fi;
     od;
     ll6:=factor(numer(add(kk,kk=ls1)));
     eqcoefs:=[op(eqcoefs),ll6];
   od;
 _EnvExplicit:=true;
   ll2:=map(x->x=0,convert(eqcoefs,set));
   ll3:={seq(op(0,op(kk,ls0))[0](op(lr0)),kk=1..ll0)};
   ll4:=[solve(ll2,ll3)];
   for jj from 1 to ll0 do
     ll2:=op(0,op(jj,ls0));
     ll4:=remove(has,ll4,ll2[0](op(lr0))=0);
   od;
   ll5:=op(ii,lr2);
   for jj from 1 to nops(ll4) do
     ll6:=op(jj,ll4);
     tm0:=[op(tm0),ll6];
     tm1:=[op(tm1),ll1];
     tm2:=[op(tm2),ll5];
   od;
  od;
 `wkptest/ltm0`:=tm0;
 `wkptest/lord`:=tm1:
 `wkptest/lpw`:=tm2:
 if nops(tm0)>0 then
    for ii from 1 to nops(tm0) do
      printf(`The leading analysis of %a branch gives:\n`,ii);
      print(op(op(ii,`wkptest/lord`)));
      print(op(op(ii,tm0)));
   od;
 else
   `wkptest/passflag`:=0;
 fi;
 printf(`\n`);
 end:

 trunform:=proc()
 global `wkptest/eqnum`,`wkptest/func`,`wkptest/var1`,`wkptest/lord`,
 `wkptest/ltm0`,`wkptest/lpw`,`wkptest/trunls`,`wkptest/lv`;
 local ii,jj,kk,nn,tm0,tm1,tm2,tm3,ls0,ls1,ls2,ll0,ll1,ll2,
 ll3,ll4,lr0,lr1,lr2,lr3,lr4,lr5,lr6,rulef,rulef1,rulef2,
 eqtruns,coef_jj,trlsii,pflag;
 `wkptest/trunls`:=[];
 ll0:=`wkptest/eqnum`;
 lr0:=op(`wkptest/var1`);
 ls0:=`wkptest/func`;
 tm1:=`wkptest/ltm0`;
 tm2:=`wkptest/lord`;
 tm3:=`wkptest/lpw`;
 for ii from 1 to nops(tm1) do
   ll1:=convert(op(ii,tm1),list);
   ll2:=op(ii,tm3);
   ll3:=op(ii,tm2);
   ll4:=-min(seq(rhs(jj),jj=ll3))-1;
   rulef1:=[];
   for jj from 1 to ll0 do
     lr1:=rhs(op(jj,ll3));
     lr2:=op(jj,ls0);
     lr3:=lr2=sum(op(0,lr2)[pp](lr0)*phi(lr0)^(pp+lr1),pp=0..ll4);
     rulef1:=[op(rulef1),lr3];
   od;
   eqtruns:=[];
   for jj from 1 to ll0 do
     lr1:=expand(eval(subs(rulef1,expand(`wkptest/lv`[jj]))));
     lr2:=combine(expand(subs(ll1,lr1)));
     eqtruns:=[op(eqtruns),lr2];
   od;
   rulef2:=convert(ll1,set);
   for jj from 1 to ll4 do
     coef_jj:={NULL};
     for nn from 1 to ll0 do
       lr1:=op(nn,ll2)+1;
       lr2:=eval(subs(rulef2,op(nn,eqtruns)))*phi(lr0)^(-lr1+1);
       lr3:=combine(expand(lr2),power);
       ls1:=[];
       for kk from 1 to nops(lr3) do
         ls2:=op(kk,lr3);
         if member(phi(lr0),{op(ls2)}) then
            ls1:=[op(ls1),ls2];
         fi;
       od;
       lr4:=factor(numer(add(kk,kk=ls1)));
       coef_jj:={op(coef_jj),lr4};
      od;
      lr1:=map(x->x=0,coef_jj);
      lr2:={seq(op(0,op(kk,ls0))[jj](lr0),kk=1..ll0)};
      _EnvExplicit:=true;
      lr3:={solve(lr1,lr2)};
      pflag:=1;
      if nops(lr3)<> 0 then
         lr4:=convert(op(1,lr3),list);
         rulef2:=rulef2 union op(1,lr3);
      else
         printf(`This branch cannot be truncated at constant term!\n`);
         pflag:=0;
         break;
      fi;
    od;
    if pflag=0 then
       next;
    fi;
    trlsii:=[];
    for jj from 1 to ll0 do
      lr1:=rhs(op(jj,ll3));
      lr2:=op(jj,ls0);
      lr3:=op(0,op(jj,ls0));
      lr4:=sum(lr3[qq](lr0)*phi(lr0)^(qq+lr1),qq=0..-lr1-1);
      lr5:=eval(subs(rulef2,lr4));
      lr6:=lr2=lr5+lr3[-lr1](lr0);
      trlsii:=[op(trlsii),lr6];
    od;
    printf(`The %a th truncated expansion form is:\n`,ii);
    for jj from 1 to ll0 do
      print(op(jj,trlsii));
    od;
    `wkptest/trunls`:=[op(`wkptest/trunls`),trlsii]:
 od;
 printf(`\n`);
 end:

 findres:=proc()
 global `wkptest/eqnum`,`wkptest/func`,`wkptest/var1`,`wkptest/lord`,
 `wkptest/lpw`,`wkptest/ltm0`,`wkptest/renlist`,`wkptest/passflag`,
 `wkptest/lv`;
 local ii,jj,kk,ll0,ld0,ll1,ll2,ll3,ll4,ll5,ll6,ll7,lr1,lr2,lr3,lr4,lr5,lr6,tm0,tm1,tm2,rulef1,rulef2,ls0,ls1,ls2,mclead;
 `wkptest/renlist`:=[];
 ll0:=`wkptest/eqnum`;
 ls0:=`wkptest/func`;
 ld0:=`wkptest/var1`;
 ll1:=`wkptest/lord`;
 ll2:=`wkptest/lpw`;
 ll3:=`wkptest/ltm0`;
 rulef1:=diffrule(ld0);
 tm0:=[];
 tm1:=[];
 tm2:=[];
 for ii from 1 to nops(ll3) do
   ll4:=op(ii,ll3);
   ll5:=op(ii,ll1);
   ll6:=op(ii,ll2);
   rulef2:=[];
   for jj from 1 to ll0 do
     ll7:=rhs(op(jj,ll5));
     lr1:=op(ld0);
     lr2:=op(jj,ls0);
     lr3:=lr2=op(0,lr2)[0](lr1)*phi(lr1)^(ll7)+lmr[jj](lr1)*phi(lr1)^(rr+ll7);
     rulef2:=[op(rulef2),lr3];
   od;
   mclead:=Matrix(ll0,ll0);
   for jj from 1 to ll0 do
     lr1:=op(jj,ll6);
     lr2:=expand(`wkptest/lv`[jj]);
     lr3:=expand(eval(subs(rulef2,lr2)));
     lr4:=expand(eval(subs(ll4,lr3)));
     lr5:=combine(eval(expand(subs(rulef1,lr4))),power);
     ls1:=[];
     for kk from 1 to nops(lr5) do
       ls2:=op(kk,lr5);
       if member(phi(op(ld0))^(rr+lr1),{op(ls2)}) then
          ls1:=[op(ls1),ls2];
       fi;
     od;
     lr6:=factor(numer(add(kk,kk=ls1)));
     for kk from 1 to ll0 do
       mclead[jj,kk]:=factor(coeff(lr6,lmr[kk](op(ld0))));
     od;
   od;
   lr1:=LinearAlgebra[Determinant](mclead);
   if lr1=0 then
      `wkptest/passflag`:=0;
       printf(`The %a th branch fails in determining the values of resonances\n`,ii);
      next;
   fi;
   lr2:=[solve(factor(lr1)=0,rr)];
   if nops(lr2)=0 then
      `wkptest/passflag`:=0;
      printf(`The %a th branch fails in determining resonances!\n`,ii);
      next;
   else
      printf(`The resonances of %a th branch are\n`,ii);
      print(op(lr2));
      printf(`\n`);
     if algbranchk(lr2)=1 then
        `wkptest/passflag`:=0;
        printf(`This branch possesses non-integer resonances!\n`);
        next;
     fi;
   fi;
   ##2003-09-03
   ##if pbranchk(lr2)=1 then
     ## `wkptest/renlist`:=[lr2,op(`wkptest/renlist`)];
     ## tm1:=[ll5,op(tm1)];
     ## tm2:=[ll6,op(tm2)];
     ## tm0:=[ll4,op(tm0)];
   ##else
      `wkptest/renlist`:=[op(`wkptest/renlist`),lr2];
      tm1:=[op(tm1),ll5];
      tm2:=[op(tm2),ll6];
      tm0:=[op(tm0),ll4];
   ##fi;
   od;
 `wkptest/ltm0`:=tm0;
 `wkptest/lord`:=tm1;
 `wkptest/lpw`:=tm2:
 printf(`\n`);
 end:

 simpltm0:=proc(tm1::set)
 global `wkptest/var1`;
 local ii,ll1,ll2,ll3,ll4,lr1,lr2,lr3,ls1,ls2,ls3,
 ls4,ls5,rule1;
 ll1:=tm1;
 ll2:=`wkptest/var1`;
 ll3:=diffrule(ll2);
 ll4:=op(remove(has,ll2,op(1,ll2)));
 lr1:={};
 lr2:={};
 rule1:={};
 for ii from 1 to nops(ll1) do
   ls1:=lhs(op(ii,ll1));
   ls2:=rhs(op(ii,ll1));
   if ls1=ls2 then
      ls3:=op(0,ls1)(ll4)=op(0,ls2)(ll4);
      lr1:= lr1 union {ls3};
      rule1:=rule1 union {op(0,ls1)(op(ll2))=op(0,ls2)(ll4)};
   else
      lr2:=lr2 union {op(ii,ll1)};
   fi;
 od;
 for ii from 1 to nops(lr2) do
   ls1:=lhs(op(ii,lr2));
   ls2:=rhs(op(ii,lr2));
   ls3:=simplify(eval(subs(ll3,ls2)));
   ls4:=simplify(eval(subs(rule1,ls3)));
   ls5:={op(0,ls1)(ll4)=ls4};
   lr1:= lr1 union ls5;
 od;
 return lr1;
 end:

 createqs:=proc()
 global `wkptest/eqnum`,`wkptest/eq`,`wkptest/var1`,
 `wkptest/func`, `wkptest/renlist`,`wkptest/lord`,
 `wkptest/ltm0`,`wkptest/lpw`,`wkptest/laureqs`;
 local tm1,tm2,tm3,tm4,tm5,tm6,tm7,ll0,ll1,ll2,ll3,ll4,
 ll5,ls0,ls1,ls2,ls3,ls4,ls5,ls6,rulef1,rulef2,ii,kk,
 jj,nn,eqs_ii,ts1,ts2,ts3,ts4,ts5,ts6,eqsr1r2;
 tm1:=`wkptest/eq`;
 ll0:=`wkptest/eqnum`;
 tm2:=`wkptest/var1`;
 tm3:=`wkptest/func`;
 tm4:=`wkptest/renlist`;
 tm5:=`wkptest/lord`;
 tm6:=`wkptest/ltm0`;
 tm7:=`wkptest/lpw`;
 ll1:=op(remove(has,tm2,op(1,tm2)));
 rulef1:=diffrule(tm2);
 `wkptest/laureqs`:=[];
 for ii from 1 to nops(tm6) do
   ll2:=op(ii,tm4);
   ll3:=op(ii,tm5);
   ll4:=simpltm0(op(ii,tm6));
   ll5:=op(ii,tm7);
   ls0:=max(op(ll2));
   rulef2:=[];
   for jj from 1 to ll0 do
     ls1:=rhs(op(jj,ll3));
     ls2:=op(jj,tm3);
     ls3:=ls2=sum(op(0,ls2)[qq](ll1)*phi(op(tm2))^(qq+ls1),qq=0..ls0);
     rulef2:=[op(rulef2),ls3];
   od;
   eqsr1r2:=[];
   for jj from 1 to ll0 do
     ls1:=tm1[jj];
     ls2:=expand(eval(subs(rulef2,ls1)));
     ls3:=-op(jj,ll5);
     ls4:=expand(eval(subs(rulef1,ls2))*phi(op(tm2))^(ls3));
     ls5:=combine(expand(subs(ll4,ls4)),power);
     ls6:=collect(ls5,phi(op(tm2)));
     eqsr1r2:=[op(eqsr1r2),ls6];
   od;
   eqs_ii:=[];
   for jj from 1 to ls0 do
     ts1:=[];
     for kk from 1 to ll0 do
       ts2:=op(kk,eqsr1r2);
       ts3:=[];
       for nn from 1 to nops(ts2) do
         ts4:=op(nn,ts2);
         if member(phi(op(tm2))^(jj),{op(ts4)}) then
            ts3:=[op(ts3),ts4];
         fi;
       od;
       ts5:=factor(numer(add(qq,qq=ts3)));
       ts6:=remove(has,ts5,phi(op(tm2)));
       ts1:=[op(ts1),ts6];
     od;
     eqs_ii:=[op(eqs_ii),ts1];
   od;
 `wkptest/laureqs`:=[op(`wkptest/laureqs`),eqs_ii];
 od;
 end:

 checkcond:=proc()
 global `wkptest/eqnum`,`wkptest/renlist`,`wkptest/func`,
 `wkptest/laureqs`,`wkptest/var1`,`wkptest/ltm0`,`wkptest/passflag`;
 local k1,ii,jj,k3,k4,k5,k6,ll0,tm0,tm1,tm2,tm3,tm4,tm5,
 ld0,eqt1,eqt2,eqt3,eqt4,eqt5,eqt6,coef_rulelist,ls1,ls2,
 rule_diff,ls0,coef_eqs,cond_k1,reson_no;
 ls1:=`wkptest/renlist`;
 ll0:=`wkptest/eqnum`;
 ld0:=`wkptest/func`;
 createqs();
 if nops(`wkptest/laureqs`)=0 then
    `wkptest/passflag`:=0;
    return;
 fi;
 ls0:=remove(has,`wkptest/var1`,op(1,`wkptest/var1`));
 for k1 from 1 to nops(`wkptest/laureqs`) do
   ##2003-09-03
   ##if `wkptest/passflag`=0 then
     ## break;
   ##fi;
   tm0:=[convert(simpltm0(op(k1,`wkptest/ltm0`)),list)];
   cond_k1:={NULL};
   coef_rulelist:={seq(op(tii),tii=tm0)};
   eqt1:=op(k1,`wkptest/laureqs`);
   for jj from 1 to nops(eqt1) do
     tm2:=op(jj,eqt1);
     eqt3:={NULL};
     for k3 from 1 to ll0 do
       tm4:=eval(subs(coef_rulelist,op(k3,tm2)));
       if tm4<>0 then
          eqt3:={op(eqt3),numer(tm4)};
       fi;
     od;
     if nops(eqt3)=0 then
        tm5:=seq(op(0,op(qq,ld0))[jj](op(ls0))=op(0,op(qq,ld0))[jj](op(ls0)),qq=1..ll0);                                            coef_rulelist:={op(coef_rulelist),tm5};
        tm0:=[op(tm0),[tm5]];
        next;
     fi;
     eqt4:=map(x->x=0,eqt3);
     _EnvExplicit:=true;
     eqt5:={solve(eqt4,{seq(op(0,op(qq,ld0))[jj](op(ls0)),qq=1..ll0)})};
     if nops(eqt5)>=1 then
       eqt6:=op(1,eqt5);
       coef_rulelist:={op(coef_rulelist),op(eqt6)};
       tm0:=[op(tm0),convert(eqt6,list)];
     else
       `wkptest/passflag`:=0;
       reson_no:=jj;
       cond_k1:={NULL};
       if nops(eqt4)=1 then
          cond_k1:=eqt4;
       else
          cond_k1:=reducecond(eqt4,jj);
       fi;
     fi;
     if `wkptest/passflag`=0 then
       break;
     fi;
   od:
 printf(`The coefficients of the %a th branch are\n`,k1);
 for k6 from 1 to nops(tm0) do
   print(op(op(k6,tm0)));
   printf(`\n`);
 od;
 if nops(cond_k1)>0 then
    printf(`The compatibility condition(s) does(do) not hold at resonance j=%a is(are):\n`,reson_no);
    print(op(cond_k1));
 fi;
 od;
 end:

reducecond:=proc(eql1::set,xx::integer)
 global `wkptest/func`,`wkptest/var1`;
 local ii,jj,ll0,ld0,lw0,lw1,lw2,lw3,lw4,lw5,lw6,
 lw7,lw8,lw9,lw10,tm1,conds,eqfs;
 ll0:=`wkptest/func`;
 ld0:=`wkptest/var1`;
 lw1:=eql1;
 lw2:=xx;
 lw3:=op(remove(has,ld0,op(1,ld0)));
 lw0:=seq(op(0,op(qq,ll0))[lw2](lw3),qq=1..nops(ll0));
 conds:={};
 eqfs:={};
 for ii from 1 to nops(lw1) do
   lw4:=op(ii,lw1);
   lw5:=indets(lhs(lw4));
   if nops(lw5 intersect {lw0})=0 then
     conds:=conds union {lw4};
   else
     eqfs:=eqfs union {lw4};
   fi;
 od;
 if nops(eqfs)<=1 then
   return conds;
 fi;
 for ii from 1 to nops(eqfs) do
   lw4:=op(ii,eqfs);
   lw5:=indets(lhs(lw4));
   tm1:=0;
   for jj from 2 to nops(eqfs) do
     lw6:=op(jj,eqfs);
     lw7:=indets(lhs(lw6));
     lw8:=lw5 intersect lw7  intersect {lw0};
     if nops(lw8)=0 then
       conds:=conds union {lhs(lw4)=0,lhs(lw6)=0};
     else
       lw9:=op(1,lw8);
       lw10:=eval(subs(solve(lw4,{lw9}),lhs(lw6)));
       if lw10<>0 then
          conds:=conds union {lw10=0};
      fi;
     fi;
    od;
 od;
 return conds;
 end:
end module:
savelib(wkptest);
